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Abstract 



Recent developments in molecular theories and simulation of ions and po- 
lar molecules in water are reviewed. The hydration of imidazole and imida- 
zolium solutes is used to exemplify the theoretical issues. The treatment of 
long-ranged electrostatic interactions in simulations is discussed extensively. 
It is argued that the Ewald approach is an easy way to get correct hydra- 
tion free energies in the thermodynamic limit from molecular calculations; 
and that molecular simulations with Ewald interactions and periodic bound- 
ary conditions can also be more efficient than many common alternatives. 
The Ewald treatment permits a conclusive extrapolation to infinite system 
size. Accurate results for well-defined models have permitted careful testing 
of simple theories of electrostatic hydration free energies, such as dielectric 
continuum models. The picture that emerges from such testing is that the 
most prominent failings of the simplest theories are associated with solvent 
proton conformations that lead to non-gaussian fluctuations of electrostatic 
potentials. Thus, the most favorable cases for second-order perturbation the- 
ories are monoatomic positive ions. For polar and anionic solutes, continuum 
or gaussian theories are less accurate. The appreciation of the specific defi- 
ciencies of those simple models have led to new concepts, multistatc gaussian 
and quasi-chemical theories, that address the cases for which the simpler the- 
ories fail. It is argued that, relative to direct dielectric continuum treatments, 
the quasi-chemical theories provide a better theoretical organization for the 
computational study of the electronic structure of solution species. 
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1. Introduction 

Water, the most commonly encountered liquid, exerts both chemical and physical influ- 
ences on aqueous molecular processes. Hydration effects are often divided into hydrophobic 
and hydrophilic categories. Hydrophilic solutes are typically ionic or polar species and may 
participate in chemical interactions with the water solvent. Because of the long range of the 
electrostatic interactions and their strength relative to k^T, hydrophilic hydration presents 
distinctive conceptual and practical issues for understanding and predicting the influence of 
hydration on chemical and biochemical events in water. 

A principal and long-standing technical issue is the treatment of infinitely long-ranged 
interactions in the context of a sample of finite size.S Recent work has helped to resolve this 
problem. One algorithmic approach to treatment of long-ranged interactions is the use of 
Ewald interactions within the conventional periodic boundary conditions.! We argue here 
that the Ewald approach is an easy way to get correct hydration free energies from molecular 
calculations, that is, to achieve well characterized results appropriate to the thermodynamic 
limit in which the system size tends to infinity for given densities and temperature. What 
is more: molecular simulations with Ewald interactions and periodic boundary conditions 
can also be more efficient than rougher approximations that are often employed to compute 
hydration free energies for molecularly well-defined problems. We anticipate results below by 
noting that we obtain accurate, thermodynamic limiting results for the hydration free energy 
of imidazole with as few as 16 water molecules included in the simulation. The price to be 
paid for this accuracy and efficiency is additional effort in understanding Ewald calculations 
from a physical viewpoint and in implementing Ewald interactions,! its equivalentsjlhi and 
alternatives.li"0 

The physical issues motivating simulation calculations of this type revolve around dielec- 
tric continuum models of hydration of ionic and polar solutes.0 It is natural and common for 
a simplified approximation to provide a conceptual baseline for considering more accurate 
theoretical results. But the converse comparison is foremost for this work. The theoretical 
efforts over recent years have provided sharper tests of the validity of the continuum ap- 
proach than merely: is an empirically correct hydration free energy obtained? Recent work 
has clarified that the dielectric models are simple implementations of thermodynamic per- 
turbation theory through second-order in electrostatic coupling parameters such as solute 
charges;0Hli dielectric models can also be considered a simple implementation of an ansatz 
that electrostatic potential fluctuations are distributed according to a eaussian probability 
density,^! or they can be considered a simplified linear response theory.&^i 

Second-order perturbation theory was found to be satisfactory for some solutes such as 
alkali ionsjl^ but unsatisfactory for waterSEl and anions.lli In the latter cases, of course, an a 
posteriori adjustment of cavity radii could still produce the correct hydration free energies.ii 
However, the more ambitious molecular theory ties the values of radii parameters to molecu- 
lar properties that depend on the thermodynamic state of the system (temperature, pressure, 
and composition of the solvent) and to non-electrostatic characteristics of the solute-solvent 
interactions. The radii are not separately adjustable when viewed from that deeper level 
of molecular theory. However, the radii can be well-defined and are not properties of the 
solutes alone but incorporate information about the solvent and thermodynamic state. 

For water as a solvent, the case of exclusive concern here, the most prominent failings of 
second-order perturbation theory are associated with solvent proton conformations that lead 
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to non-gaussian fluctuations of electrostatic potentials.oEl Thus, the most favorable cases 
for the second-order perturbation theories are classic positive ions. In such cases, oxygen- 
hydrogen bonds of water are oriented away from the ion. Neutral, polar molecules that may 
form specific hydrogen bonds with the solvent are more challenging for these theories, though 
the hydration free energies sought are smaller in magnitude than for typical ions. Negative 
molecular ions are expected to offer further complications because now the problematic 
proton interactions with the solute will be strong. However, we have less experience with 
realistic negative ions partly because the molecular models used for simulation are less well 
developed than for other cases. 

The appreciation of the different possibilities for fluctuations has led to new theories 
of electrostatic hydration free energies.iil These theories analyze electrostatic distributions 
more broadly, still using gaussian models at crucial steps; but now several gaussian distri- 
butions are derived from an analysis of the first shell environment of the solute. For the 
important case of hydration of a water molecule, this extension repairs the breakdown of 
a single gaussian theory. Negative ions can still be problematic but the multiple gaussian 
approach has also motivated development of quasi-chemical theorieJ^ that are based, in 
principle, on full information about the thermal motion of the first hydration shell. Though 
experience with the quasi-chemical theories is limited,ii we anticipate that they should pro- 
vide better descriptions of the hydration free energies, in addition to providing a reasonable 
pathway to carry-out solution phase electronic structure calculations on hydrated negative 
ions, calculations that would be difficult particularly in the absence of hydration effects. 

In the following section, we will first introduce the model solute imidazole, which was cho- 
sen as a molecular solute to exemplify, combine, and extend aspects of ionic and polar solutes 
studied previously.llHh0'0@ Results for imidazole and imidazolium will be used throughout 
the manuscript to illustrate the theoretical issues. We will discuss the Ewald treatment of 
electrostatic interactions, motivating it in various ways. Subsequently, finite-size effects will 
be studied. The correction for the typically large finite-size effects is essential for accurate 
calculations of solvation free energies of polar and charged solutes. We distinguish between 
electrostatic finite-size effects that are independent of the thermodynamic state and the 
characteristics of the solute, and the remaining thermodynamic finite-size effects. We will 
then introduce perturbative methods for calculating solvation free energies that are based 
on the approximately gaussian character of the electrostatic potential fluctuations. Non- 
gaussian behavior and its accurate treatment using multistate gaussian and quasi-chemical 
models will be the focus of the last section. 



2. Example: Imidazole and Imidazolium in Water 

To illustrate the various issues arising in calculations of solvation free energies of charged 
and polar molecules, we present new calculations of the hydration of imidazole and imida- 
zolium. We choose this example because recent interest in these problems has focused on 
predicting acid-base equilibria of biochemical relevance.iH^ii We will calculate the charg- 
ing free energies of the protonated imidazolium and the neutral, polar imidazole (Figure p. 
Imidazole in water provides a rich example: the polar imidazole molecule can be protonated 
at the N3 position to form a molecular cation, imidazolium. This protonation reaction 
has a pK^i of about 7.0 It provides a basis for pKa calculations of ionizable residues of 
proteins.EM^ Imidazole is the building block of histidine, one of the most active amino 
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acids enzymatically and ubiouitous in the active sites of enzymes that operate at room 
temperature and neutral pH.E2l 

The protonation of imidazole has been studied previously using combinations of dielec- 
tric, quantum mechanical, and computer simulation methods. EjH Here, we will focus on 
the solvation contribution to the protonation equilibrium and reserve quantum mechanical 
intramolecular effects for subsequent treatment. 

The imidazole/imidazolium system is more complex than the systems we have studied 
before using explicit solvent models: mono- and divalent ions,0ill water in water,lli and 
tetramethylammonium.B The analysis below should also illustrate how the calculation of 
solvation free energies using the Ewald method and equilibrium fluctuations of electrostatic 
potentials can be extended to proteins, in particular the calculation of pica's of amino acids. 

Monte Carlo computer simulations of imidazole in water. We studied the sol- 
vation of imidazole [Im(p)] and imidazolium [Im(+)] in water using Monte Carlo (MC) 
simulations in the canonical ensemble. For the Im(p) and Im(+) molecules, we used the 
partial charges and geometry of Topol et al.,0 as compiled in Table |. For water, we used 
the SPC/E modeLej The temperature was 298 K. The Ewald method was used for the 
long-range electrostatic interactions with a real-space screening factor of 77 = 5.6/L, where 
L is the length of the periodically replicated cubic box. A cutoff of k'^ < 38(27r/L)^ was 
applied in Fourier space, resulting in 2 x 510 k vectors being considered. A cutoff of L/2 
was applied to the Lennard- Jones and real-space electrostatic interactions based on atoms. 
The background dielectric constant in the Ewald method was corrected from infinity to 80.0 
The partial molar volume of the imidazole was chosen as two times that of bulk water at a 
density of 997.07 kg/m'^, such that the pressuree^ was about one atmosphere. 

The Metropolis MC method was used to sample configurational space in the 
simulations.^ The translational and rotational move widths of water were chosen to give 
about 40 % acceptance ratios. The solute was allowed to move as well. Simulations started 
from random configurations or configurations of previous runs with different charges, equi- 
librated for at least 100,000 and 50,000 MC passes, respectively, where one pass is one 
attempted move for each of the particles. Electrostatic potentials at solute atom positions 
and binding energies of the solute were calculated after every pass using the Ewald method 
in simulations extending over 200,000 MC passes each. Simulations were performed for 
the uncharged, half-charged, and fully charged Im(p) and Im(+) in their respective geome- 
tries. The equilibrium simulations were performed with 16, 32, 64, 128, 256, and 512 water 
molecules to study finite size effects. 

To complete the thermodynamic cycle and check for consistency, two runs of slow growth 
thermodynamic integration were used to calculate the free energy of converting the geometry 
from the uncharged Im(p) to the uncharged Im(+) conformation within 150,000 MC passes. 
Six runs of 200,000 MC passes were used to calculate the free energy of converting the polar 
Im(p) into an Im(-(-) cation, starting from different equilibrated configurations and averaging 
three charging and three uncharging runs. The thermodynamic integrations were carried 
out with 256 water molecules. 

We will discuss the results of these calculations as they come-up in the theoretical nar- 
rative. However, before considering more subtle issues we can make a direct comparison 
of the average electrostatic potential exerted by the solvent observed during the simulation 
with the corresponding predictions of dielectric models. Figure shows that comparison for 
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several sets of radii in current use. Such a comparison illustrates the basic issue of sensitivity 
of thermodynamic results to the radii parameters and whether extant empirically adjusted 
radii are transferable to slightly different cases. 

3. Noteworthy Aspects of the Ewald Treatment of Electrostatic Interactions 

Viewed from an historical perspective, the most appropriate treatment of electrostatic 
interactions for simulation calculations has been a contentious issue. In this setting it is 
helpful to note some broad, and non-technical, characteristics of Ewald treatments that 
might be typically overlooked. We preface these observations by noting that simulation 
calculations treat finite systems. Most commonly, periodic (or Born-von Karman) boundary 
conditionsSS are utilized for the exterior boundary of the finite system considered. The 
theoretical issues engendered by these boundary conditions with finite-rang and long- 

ranged interactions are reasonably well understood. Of course, simulation calculations need 
not address issues of what is happening outside the simulation cell. We note that it is possible 
to compute the Ewald potential - and in more than one way - without consideration of image 
charges outside the simulation cell. So intuitive arguments based upon image charges can 
be avoided completely. 

It is convenient to express the Ewald electrostatic energy of a system of partial charges 
Qia at positions rjQ, on molecules z as a sum of effective pair interactions and self terms, 

i,j a,f3 



J2 J2 



if3 



I 1 



i a 



(1) 



where Yiajjs = rj^—Tia. The Coulomb energy U was split into intermolecular, intramolecular 
and self interaction contributions. The Fourier representation of ip{r) reveals the periodicity 
of this potential^ 

V'(r) = ^ E P"--' ■ (2) 

The k sum extends over the reciprocal lattice derived from the real-space lattice n of pe- 
riodically replicated simulation boxes. For a cubic lattice of length L = V^^^, we have 
n = L{i,j,k), and k = 27rL~^ {i,j,k), where i, j, and k are integers. For numerical con- 
venience, ip{r) is partly transformed into real space, which leads to its Ewald lattice sum 
representation: 

_ ^ erfc(^7|r + n|) , >^ 47r .2/4^2+,^^ ^ /o^ 
^^'^"il^ |r + n| +4^V^" W 

1] is a. convergence parameter that is chosen to accelerate numerical convergence. Note that 
the value of ip{r) is independent of 
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and that the average potential in the box is zero,0'iiliiiil 

(^(k = 0)= [ dvip{r) = . (5) 

In the following, we will separate the Coulomb energy Ud of the solute from the total 
Coulomb energy eq When the system contains a solute with partial charges qp at positions 
r/3, its electrostatic interaction energy can be split into a solvent term and self- interactions, 

^el = 51 QlSQiay^i^fSia) 



+ J2 mi 
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+ - qs^ lim 



1 



(6) 



The first, second, and third sum are the direct interactions with water, the interactions of 
charges on the solute with other solute charges, and the self-interactions of solute charges, 
respectively. 

The Ewald potential ^^(r) is the solution of the Poisson equation that is peri- 
odic with the fundamental period of the simulation cell.llli3ilii A periodic solution 
of the Poisson equation requires that the surface integral of the electric field normal to the 
surface of the simulation cell be zero. This means that the material in the simulation cell 
must have zero net electric charge. If the physical system of interest is non-neutral, a uni- 
form background distribution of charge is included to neutralize the non-zero charge of the 
physical system. 

Consider an elementary charge. Because the Ewald potential is periodic, we can consider 
the Ewald electrostatic potential implied by centering the simulation cell on this elementary 
charge. By symmetry the Ewald normal electric field is zero on the cell boundary. The 
Ewald potential can thus be considered to be that of a cut-off on the cell boundary - the 
cut-off at the maximum distances achievable - with zero normal derivative analogous to a 
shifted-force correction. 

The Ewald potential pushes the electrostatic boundary outward as far as 
possible but still retains smoothness on the boundary. The minimum image cut-off 
shares with the Ewald treatment the property that the electrostatic potential is not cut-off in 
any region of the simulation volume, the largest volume that must be physically considered. 
However, as demonstrated above, the Ewald potential is smooth on that boundary since it is 
the periodic solution of the Poisson equation. This is an important technical advantage that 
facilitates investigation of system size dependence of computed properties, i.e., the variation 
of system properties with variations of the cell boundaries. In those cases where the physical 
system of interest is non-neutral it is a helpful point of view that the background charge 
density is a simple device that permits smoothness of the computed potential on the system 
boundary. We emphasize that the effect of the neutralizing background charge disappears as 
the thermodynamic limit is approached, in which the background charge density disappears. 
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Reaction field potentials are Ewald potentials for different background 
charges. Reaction field methodsi'iiUlSl are computationally efficient alternatives to Ewald 
summation. The effective potentials of the site-site reaction field (SSRF) methodH and the 
generalized reaction field (GRF) method0 can be viewed as Ewald potentials for a non- 
homogeneous background charge. The SSRF and GRF potentials are the solutions to the 
Poisson equation for a source charge and a compensating background. The background 
charge densities in the SSRF and GRF method are that of a homogeneous sphere and a 
radially symmetric charge distribution centered around the source charge, respectively. The 
SSRF and GRF method have a finite range because of the radially symmetric background 
charge densities that exactly compensate the source charge when the cut-off distance is 
reached. 

Comparison of Ewald potentials from simulation of water to electrostatic 
potentials in isolated water droplets. It is interesting to make some simple numerical 
comparisons between the Ewald potentials that are experienced in simulation of water with 
periodic boundary conditions to the corresponding electrostatic potentials in water droplets. 
Such a comparison is simplified if we locate a distinguished solute at the origin of our 
Cartesian coordinate system. For a spherical Lennard- Jones solute, Figure ^ shows a typical 
variation of the electrostatic potential at the solute center with inclusion of the charge 
density in progressively larger spherical volumes of radius R around the solute. Notice the 
substantial variation of the electrostatic potential with inclusion of the solvation shells near 
the solute. However, after about three shells the net electrostatic potential oscillates about 
the Ewald asymptotic value before the ball penetrates the physical interface of the droplet. 
Thereafter, the net electrostatic potential displays the effects of the surface polarization of 
the droplet as it makes a transition to the very different value that characterizes the whole 
droplet. It is clear in this case that the Ewald potential faithfully captures this interior 
potential while avoiding detailed considerations of the droplet interface. 

Single ion hydration free energies are well-defined within molecular simu- 
lations. Electrostatic potentials can be defined unambiguously as solutions of the Pois- 
son equation with specified charge densities and boundary conditions. Electrostatic poten- 
tials are computed throughout simulations of aqueous solutions even if charge densities and 
boundary conditions may not be specified explicitly. If the system of interest is non-neutral, 
these issues deserve emphasis because single ion free energies are typically not measured 
experimentally. 

On a molecular scale, dependences on specifics of the boundary conditions in the defini- 
tion of single ion hydration free energies can be avoided. This is accomplished by spherically 
integrating the electrostatic potential over the charge density around a charge site up to 
a distance where thatnotential saturates. Using Ewald interactions corresponds to such a 
spherical integration.cllEj Figure ^ shows results for electrostatic potentials obtained using 
two different choices for the boundary conditions - a solute-water cluster and a periodic 
system. The observed agreement is a non-trivial computer experimental observation. 

4. System Size Extrapolation 

Computer simulations are performed for a finite system of molecules. In most applica- 
tions, the properties of the thermodynamic-limit, infinite system are sought. In coulombic 
systems, pronounced finite-size effects are ubiquitous due to the long range of the interac- 
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tions. We can separate finite-size effects in coulombic systems into two categories: (1) those 
caused purely by the long-range electrostatics independent of the thermodynamic state; and 
(2) finite size effects that depend on the thermodynamic state (temperature, pressure, etc.). 

The electrostatic finite-size effects in the Ewald treatment of Coulomb interactions arise 
from self-interactions and interactions with the neutralizing background essential for non- 
neutral systems under periodic boundary conditions. Electrostatic finite-size effects can be 
treated exactly by including the second and third sum in the electrostatic energy Uei of the 
solute eq ^, which account for self-interactions of the solute.cJO For an ion of charge g, 
the resulting correction to the solvation chemical potential is 

f^elec A'-sim 2 ' 

where /Xsim is the chemical potential for charging the ion from zero charge to net charge 
q calculated from the Ewald interactions with the solvent excluding self-interactions {i.e., 
including only the first sum in eq fieiec includes the self-interactions; and ^ is the ionic 
self-term. For a cubic box of length L, we have ^ = limr^o[v'(i") ~ Vkl ~ — 2.837297/L. 
Electrostatic finite size corrections for polar molecules are developed in Ref. ^ The corre- 
sponding free energy of changing partial charges located at positions r„ on a molecule from 
g« to q'^ is 



^ E - e (8) 



2 

where (. . .) denotes a canonical average. A/isim includes only the energy difference corre- 
sponding to the first sum in eq ^ excluding self-interactions. Note that in constant pressure 
simulation with varying box volume, ^ must also be averaged. 

Thermodynamic finite-size effects on the other hand can only be corrected approximately 
within a model. For instance, we can use the difference between a truly infinite version of 
a model and its finite periodic version to correct for thermodynamic finite size effects,0H^ 
as schematically shown in Figure |[ For spherical ions, a Born modelS and its periodic 
equivalent leads to finite size corrections that depend on the dielectric constant e of the 
solvent, an effective Born radius Rb of the ion, and the net charge q of the ion:Ey't3 



1 2 

/^therm ~ /^elec ~^ 2 ^ 



47r(6-l)i?B^ - 
e + 3eL3 ' 



where /xtherm is the chemical potential for charging that includes the thermodynamic and 
electrostatic finite-size corrections. 

Figure § illustrates the finite size effects for the free energies of charging Im(p) and Im(-(-) 
after correction for Ewald self-interactions. Free energies were calculated from sixth-order 
integration formulas with corrected means and variances from Table |I| that include the 
electrostatic^ but not the thermodynamic finite-size correction. The free energy is plotted 
as a function of where L is the box length. We find that the free energy of charging 
the polar Im(p) is independent of the system size within the statistical errors of about 1 
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kJ/mol for = 16 to = 512 water molecules. However, the free energy of charging the 
Im(+) cation shows a system size dependence proportional to l/I/^, as would be expected 
from our finite size analysis-S Rather than using a more realistic shape of the molecule in 
the dielectric model of ref |60|, we fit the observed free energies to to the spherical Born 
model eq|^ with an effective radius Rb = 0.207 nm and an infinite dielectric constant which 
reproduces the data over the whole range of system sizes of 16 < A^ < 512 water molecules. 

To further illustrate the power of the finite-size corrections, Figure ^ shows the proba- 
bility distributions of electrostatic interaction energies f/ei of the imidazolium cation Im(+) 
with the A^ water molecules eq ^. Also shown in Figure |^ are the corresponding gaussian 
distributions, which nicely reproduce the calculated histograms of Uei- However, we observe 
a strong system size dependence: small systems have narrower distributions of f/ei with 
less negative averages compared to large systems. When we apply the electrostatic and 
thermodynamic finite size corrections for the mean and average, the gaussian distributions 
"collapse" to a single distribution corresponding to the limit of an infinite system size. We 
note that the finite-size correction is large: for A^ = 16 water molecules, the average Uei 
changes by about —370 kJ/mol (250 k-QT). 

5. Perturbation Theory 

A fundamental view of the thermodynamics due to electrostatic interactions may be ob- 
tained from consideration of the distribution p{u\ A = 0) of electrostatic energies u in the 
reference charge state A = 0. The part of the chemical potential due to electrostatic inter- 
actions A/i(A), the thermodynamic parameter sought is then expressed by the fundamental 
result 

g-/3AMA) _ ^e-^^")^_^ = j du p{u; A = 0)6"^^". (10) 

Here, (3~^ = k-^T is Boltzmann's constant times the temperature and (. . .)a=o denotes a 
thermal average with the solute in reference state A = 0. This formula requires the consid- 
eration of the electrostatic potential even though the electrostatic potential of a phase is an 
operationally subtle property.ilii Despite that subtlety, the potential sought is conceptually 
well-defined as the solution of Poisson's equation with specified charge density and boundary 
conditions.0ilii'§ 

Direct use of eq ITUI can present difficulties. Though p{u) is often substantially gaussian. 



the fundamental formula eq|Ty is sensitive to the tails of p{u). That limits the applicability 
of eq|lO|for calculations of even small changes in the charge state A. In addition, the simple 
estimator ln(e~^'^")A=o ~ ln[M~^ X^i^i e~^^"'] from M energies Ui observed in a simulation 
is biased and large sample sizes M are required for this bias to be negligible.i3 

Perturbation or cumulant expansions provide a technique to analyze these 
distributions.EHlliUli A cumulant expansion^ with respect to A of eq |10| provides 



{exp{-p\u))^^Q =exp 



oo 

1^0 ^! . 



(11) 



This defines the cumulants C„ of order n = 0, 1, 2 as 
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Co = (12a) 
Ci = («),=o (12b) 
C2 = {{u-{u)^^,f)_^ . (12c) 



A=0 



We interpret eq |TT] as a Taylor expansion in A but augment A/i(A) to include the self- 
contribution (Ag)^^/2, where ^ vanishes in the thermodynamic limit but accounts for finite- 
size effects as discussed above. Then for the charging of an ion from a neutral reference 
condition we have 



(Ag) 



2 



2 



A/i(A) = Ag(«),^o-^ f3{{u-{u),^,y}^^^-^ +■■■ . (13) 



This result should be compared to the Bornll formula for the hydration free energy due to 
electrostatic interactions of a spherical ion of radius R and charge Ag: 

where e is the dielectric constant of the solvent. The matching of the second-order terms 
between the cumulant expansion and the continuum formula provides a determination of the 
radius R of a spherical ion. The distributions, p{u; A = 0), required to evaluate the cumulant 
averages involve the non-electrostatic interactions, A=0, between solvent molecules and the 
solute. The indicated average thus generally depends upon full characterization of the 
solvent. Note that the continuum model neglects the molecular contribution linear in A. 
This linear term contributes to the asymmetry between anion and cation solvation, making 
the solvation of anions more favorable for a given ion size.0 

In principle, higher-order cumulants could be used to obtain information about the other 
Taylor coefficients. However, as was observed by Smith and van GunsterenJ^l higher-order 
cumulants are increasingly difficult to extract from computer simulations of limited duration. 
Though direct extension of perturbation theory beyond fourth order has been impractical, 
interpolative approximations polynomial in A have been more successful. For the charging 
of water and ions,_nolynomials of order six and higher were necessary to account for the 
simulation data.li3llil'0 Thus, perturbation theory was found to be unsatisfactory in such 
cases. For atomic0 and molecular ions," a kink is typically observed for dA^(X)/dX as a 
function of charge A at modest values of this parameter when the solvation shell changes 
from a cationic to an anionic structure. Additional nonlinearities were observed at high 
values of the ionic charge A. 

Table |T| contains the averages and variances corrected for electrostatic finite-size effects 
of the electrostatic energies of Im(p) and Im(+) for system sizes between 16 and 512 water 
molecules. Errors of one standard deviation of the mean were estimated by plotting the 
block error as a function of the number of blocks. The estimated error reaches a plateau 
when block values are uncorrelated. From the averages and variances, we can calculate the 
chemical potentials of charging using integration formulas (ijk) exact to various orders that 
involve i, j, and k derivatives of the free energy with respect to the coupling parameter at 
the uncharged, half-charged, and fully charged state:0 

A/x(010) ^ Ci(A = 0.5) (15a) 
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~60 [^^(^ = 0) - ^^(^ = 1)] (15e) 

where the cumulants Ci and C2 contain the electrostatic finite-size corrections. Formulas 
involving higher cumulants and different nodes Aj are discussed in Ref. |7l|. The integration 
formulas eqs_| 15a| to [15e| are exact to order 2, 2, 4, 4, and 6 in a perturbation expansion, 
respectively.EyOFigure |^ shows the free energy difference between Im(-|-) and Im(p) as a 
function of the integration order for the 512 water molecule system. We find that as the 
order of the integration formula increases, the free energy difference converges, with the 
sixth order formula bracketed by the two fourth-order formulas. The statistical error of the 
free energy difference is about 1.5 kJ/mol. Notice that the discrepancy between the two 
second-order results of Figure ^ is significant on the scale of the statistical uncertainties. 
This emphasizes that the charging free energy is not a Quadratic function of the coupling 
parameter. Note that the centered second-order formula!!^ has a smaller systematic error. 

Figure || illustrates the complete four-node thermodynamic cycle, where A is a coupling 
parameter changing the partial charges on the molecule linearly from state zero to one. 
The four nodes of the cycle are the uncharged and charged imidazole and imidazolium. We 
find that the free energies of charging and conformational changes are consistent within the 
statistical errors. Interestingly, the free energy of charging the polar Im(p) to the Im(+) 
cation has a maximum for the linear charging path chosen here. This increase refiects the 



linear terms of eq |T3[ i.e., increasing the net charge on the imidazole initially costs free 
energy. 

Dielectric continuum models predict a quadratic proportionality of the free energy of 
charging on the linear coupling parameter A. In a molecular theory, such a quadratic charg- 
ing free energy arises when the probability density of electrostatic potential fiuctuations is 
gaussian. Second-order perturbation theory would then be exact. Figure ^ compares second- 
order perturbation theory with the reference sixth-order free energy polynomial calculated 
from the averages and variances in Table |l[ We find that the perturbation expansions about 
the charged state (A = 1) are accurate over a relatively wide range from A = 1 to almost 
A = 0.2. The expansion about the uncharged state A = on the other hand breaks down 
rapidly at about A = 0.2. 



6. Non-Gaussian Fluctuations 

Multistate Gaussian Models. One idea for improvement of dielectric models is based 
upon a physical description of the structure of the first hydration shell. It can be viewed from 
the perspective of Stillinger- Weber inherent structures or substates.0 These are potential 
energy basins of attraction for steepest-descent quenching of first hydration shell molecules. 
If those first hydration shell molecules stayed always in one basin, then a gaussian model 
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for thermal fluctuations would be reasonable. Empirical radii parameters reflect the char- 
acteristics of that single basin. However, changing conditions may result in reweighting of 
slightly accessible basins or the opening of new basins. The gaussian or dielectric models 
may fail to describe these possibilities well. This picture is physically better defined than 
the commonly nonspecific discussions of electrostriction and dielectric saturation. 

A corresponding "multistate gaussian model" was developed in ref ^ Attention is di- 
rected to the thermal probability distribution of electrostatic potential energies of the solute. 
Rather than approximating this distribution as a single gaussian distribution, perhaps with 
perturbative corrections, we discriminate hydration structure on the basis of simple pa- 
rameters diagnostic of hydration substates. We assume that the probability distribution of 
electrostatic potential energies is gaussian for each substate. Therefore the full distribution 
is a superposition of gaussian distributions for the various substates. 

Thus we attempt to represent the observed complicating features of p{u) by a combination 
of simpler states: 

= Yj'^riPniu) , (16) 

n 

with weights Wn > 0,J2n'^n = 1 and normalized densities Pn{u) > 0, J du Pn{u) = 1. We 
will seek pn{u)'s of gaussian form, representing the overall system as a linear combination of 
gaussian subsystems, each showing linear response to electrostatic interactions. Representing 
p{u) by a sum of gaussian densities can give nontrivial results for the chemical potential, as 
can be seen by substituting eq |T^ into eq |10|, 

A/i(A) = -ksT In yj^e-f'^n^„+f3-^x-^aJ/2 ^ ^^^^ 

n 

where m„ and are the mean and variance of the gaussian p„, respectively. 

The non-gaussian fluctuations of the electrostatic potential in liquid water are associated 
with changes in the conformations of protons that make hydrogen bonds to the solute. If 
those fluctuations could be tempered, a gaussian model might become more accurate. Thus, 
suitable substate diagnostic parameters are the number of hydrogen bonds made to the 
solute. 

Explicit calculations have shown that this approach eliminates most of the detailed nu- 
merical inaccuracies of the gaussian fluctuation models for hydration of a water molecule 
in liquid water.ill The markedly non-gaussian p{u) was accurately represented as the sum 
of gaussian distributions implied by this definition of a hydration substate. We found 
Wn > 10~^ for 1 < n < 6 with 3.64 being the average number of neighbors and n = 4 
the most probable number of neighbors. The calculated change of the chemical potential 
upon change of the charge state of a solute water molecule is correct to within 5 %. This is 
a remarkable result because A/i(A) is non-quadratic, requiring an eighth-order polynomial 
to fit the simulation data for chemical-potential derivatives.ll3oO This shows that sufficient 
information can be extracted from the simulation to describe the distribution p{u) helpfully; 
and that such an approach can be successful even for perturbations involving changes of the 
chemical potential as large as 14 fee 7". 

Similar behavior can be anticipated for hydration of other neutral, polar solutes such 
as the imidazole example studied here. Figure ^ (inset) shows the results of the multistate 
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gaussian model applied to charging and uncharging the polar imidazole Im(p). Fluctuation 
data were collected from a simulation of the uncharged and charged Im(p) in = 128 water 
molecules, extending over 10^ MC passes to allow for error estimates. Instead of determining 
the overall mean and variance of the electrostatic potential for a second-order perturbation 
expansion, we calculate the means and variances for several gaussian distributions from 
structures sorted according to the number of hydrogen bonds. Inspection of the radial 
distribution functions of water oxygen and hydrogen around imidazole sites shows one strong 
hydrogen bond donor, HI, and one acceptor, N3. As a criterion for the formation of a water- 
imidazole hydrogen bond, we used that the distance between the acceptor (nitrogen N3 or 
water oxygen) and the donor hydrogen has to be smaller than 0.23 nm. We can then sort 
structures according to the numbers of hydrogen bonds accepted and donated by the solute. 
With this simple criterion, we find that six and two gaussian distributions contribute to the 
expansions about the charged and uncharged state of Im(p), respectively. This multistate 
gaussian model greatly improves the quality of the expansions, both about the charged and 
uncharged state. The reference free energy is now within the statistical errors of the two 
multistate gaussian models over the whole range < A < 1. 

Quasi-chemical theories. Those more difficult anionic cases mentioned above can be 
attacked more directly. The local neighborhood is again used to discriminate structural 
possibilities. But, in addition, the consequences for the hydration free energy of the molecu- 
lar interactions within that neighborhood are treated fully. This reserves the longer-ranged 
interactions for simple approximations, e.g. with gaussian models. 

These theoretical developments arose from recent molecular calculationsil that suggested 
how a chemical perspective can be helpful in computing thermodynamic properties of water 
and aqueous solutions. That calculation used electronic structure results on the Fe(H20)6^^ 
cluster and simple, physical estimates of further solvation effects. The results were organized 
according to the pattern of a simple chemical reaction and a surprisingly accurate evaluation 
of the hydration free energy was obtained. Despite this recent motivation, the theories devel- 
oped are akin to good approximations of historical and pedagogical importance in the areas 
of cooperative phenomena and phase transitions.El In those areas, similar approximations 
are called Guggenheim, quasi-chemical, or Bethe approximations. 

These quasi-chemical theorieS are constructed by considering a geometric volume fixed 
on the solute molecule and performing a calculation analogous to the evaluation of a grand 
canonical partition function for that volume. All the possibilities for occupancy of that 
volume by solvent or other solution species must be considered eventually. The final result 
can be described by reference to simple patterns of chemical equilibria such as: 

X"- + nH^O ^ X{H20)n''~ . (18) 

The solute of interest is denoted here generically as X'^~] a star- type cluster of that solute 
with n water (W) molecules is considered as the product. For such a cluster, call it an 
M-cluster, we could calculate the equilibrium ratio Km for a dilute gas phase. Note that 
KmPw^ is dimensionless and this observation resolves standard state issues. The factors 
denoted by (exp{— /5Am})o,c, where C indicates either a water molecule or the M-cluster, 
carry information about solvation free energy of the species involved. For the species other 
than the cluster this is the familiar Widom factor. For the cluster, this factor requires 
slight additional restriction but can be verbally defined by saying that it is the average 
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of the Boltzmann factor for cluster-solution interactions over the thermal motion of the 
cluster and solution under the condition that the only interactions between these subsystems 
rigidly exclude additional solvent molecules from the cluster volume for the complex. This 
restriction enforces a constraint required to preserve simple enumerations that underlie these 
results. The theoretical structure is designed so that simple approximations such as dielectric 
models might be used for the factors {exp{—PAu})o^c- But more detailed techniques might 
be applied to the calculation of Km- Finally, we compile 

- (exp{-/?AM})o,M 

Km = Km-T, r n^, 1\ T^T " ^^^^ 

Thus, this Km is built on the pattern of the chemical equilibrium eq |l^ but without a 
'solvation factor' for the 'bare' solute. 

Now consider all possibilities for clusters. A thermodynamic implication of this informa- 
tion is: 

fxx-^- = kT ln[px.- V/qx.-] - kT ln[po E KmPw""] • (20) 

M 

Po is the probability that the clustering volume would be observed to be empty in the 
equilibrium solution; thus — kT In po is the free energy for formation of a cavity for the 
clustering volume in the solution. The sum is over all clusters with zero or more ligands. 
The product of the densities involved with each term includes a density factor for each ligand. 
This formula makes the conventional separation between the contributions of intermolecular 
interactions and the non-interaction (ideal) terms; qxi- is the partition function of the 
bare solute in the absence of interactions with any other species and pxi- is the density of 
the solute. As an example, for an atomic ion such as the chloride ion Cl~ we would put 
Qxi- = V/A^ with V the volume of the system and A the thermal deBroglie wavelength of 
the chloride ion. This formula becomes approximate when approximate models are adopted 
for pq, for Km, and for the solvation factors. Those quantities depend on definition of the 
clustering volume. But, since the physical problem is independent of those parameters, the 
theory should be insensitive to them. 

The motivation of this approach is the fact that a substantial but intricate part of the 
free energy sought is to be found in Km- The number of possibilities for ligand populations 
will be small for molecular scale clustering volumes. So a limited number of terms must be 
considered. Because the clusters will be small systems, elaborate computational methods 
can be applied to the prediction of the Km, including current electronic structure techniques. 
With the complicated chemical interactions separated for individual treatment the remaining 
hydration contributions should be simpler and the required theories better controlled. 



Equation ^ should be compared with eq One difference is that eq ^ attempts to 



provide the whole hydration free energy, not just the part due to electrostatic interactions. 
That explains the presence of po in eq |2^. Beyond that, the structures of these formulas 
are similar. The presence of more than one term in the sum of eq 120 is an expression of 



an entropy contribution associated with the possibilities for different ligand populations. 
Finally, the complete calculation of the Km includes non-gaussian statistical possibilities 
not anticipated by eq O. 
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7. Conclusions 



Recent calculations of the hydration free energy due to electrostatic interactions between 
charged and polar solutes in water have obtained high accuracy results for the simple molec- 
ular models that are the basis of most simulation calculations.SBii'iiJliisI An important 
step in securing those high accuracy results has been a careful consideration of treatment 
of long-ranged interactions. That work suggests that the Ewald method is an easy way to 
get correct hydration free energies from molecular calculations, that is, to achieve well char- 
acterized results appropriate to the thermodynamic limit in which the system size tends 
to infinity for given densities and temperature. Additionally, it suggests that molecular 
simulations with Ewald potentials and periodic boundary conditions can have efficiencies 
comparable to rougher approximations that are often employed to compute hydration free 
energies for molecularly well-defined problems.^ And furthermore, this has produced a sim- 
ple, effective, and clear understanding of how to extrapolate electrostatic hydration free 
energies to the thermodynamic limit; an accurate evaluation of the hydration free energies 
of imidazole and imidazolium can be obtained with as few as 16 water molecules included 
in the simulation. 

These high accuracy results for well-defined models permitted careful testing of simple 
theories of electrostatic hydration free energies. The simplest theories, dielectric continuum 
models, have been found to be rough despite the fact that they can always be adjusted to 
reproduce an empirical answer given a priori. Such a conclusion has surely been widely 
expected. However, the testing has led to new theories, the multistate gaussian and quasi- 
chemical theories, that should permit more revealing molecular scale calculations. The 
quasi-chemical approaches seem to provide the most natural way to utilize current electronic 
structure packages to study electronic structure issues for solution species. This should be 
particularly helpful for treatment of basic, molecular anions that are ubiquitous in aqueous 
solution chemistry. 

Physical conclusions more specifically are that the most prominent failings of the sim- 
plest theories are associated with solvent proton conformations that lead to non-gaussian 
fiuctuations of electrostatic potentials. Thus, the most favorable cases for the second-order 
perturbation theories are monoatomic positive ions. In such cases, oxygen-hydrogen bonds 
are oriented away from the ion, placing those protons as far out as possible. Neutral, polar 
molecules that may form specific hydrogen bonds with the solvent are more difficult for 
these theories, though the hydration free energies sought are smaller in magnitude. Nega- 
tive molecular ions are expected to offer further complications because now the problematic 
proton motions occur close to the solute and the hydration effects will be larger for anionic 
species. 
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FIGURES 



Figure 1. Protonation equilibrium between imidazole and imidazolium. 



Figure 2. Comparison of dielectric models (ordinate) with molecular simulations (ab- 
scissa) for the induced electrostatic potentials due to the solvent at the atom centers 
for Im(+) (upper panel) and Im(p) (lower panel). Dielectric model results were ob- 
tained for several sets of radii in current use: diamonds: Rc=0.267 nm, Rjv=0.231 nm;@ 
large circles: Rj:/(Ar)=0.1160 nm, Rj:/(c)=0.1710 nm, Rc=0.230 nm, R7v=0.150 nm;ii small 
circles: Rh=0.1172 nm, Rc=0.2096 nm, Rjv=0.1738 nm;i crosses: Rh=0.1172 nm, 
Rc=0.1635 nm, R7v=0.1738 nmB A boundary element method was used for the dielec- 
tric model calculations.0 Notice that a radii set that happens to be qualitatively satisfac- 
tory for the cation (diamonds) can be significantly less satisfactory for the slightly different 
circumstance of the neutral polar molecule. 



Figure 3. Electrostatic potential at the center of a neutral Lennard- Jones solute in SPC 
wateril from simulations of a solute at the center of a cluster with 1024 water molecules 
(dashed line) and with periodic boundary conditions and Ewald summation (solid line). 
Shown is the potential obtained by integrating the charge density around the solute up to a 
distance R using 1/r (cluster) and (f{r) (Ewald) for the Coulomb interactions. The results 
are those of ref ESl 



Figure 4. Schematic representation of the thermodynamic finite-size correction. The ther- 
modynamic finite-size correction {. . .} is the difference between an infinite Born model and a 
Born model under periodic boundary conditions. A spherical ion of charge q and radius Rb 
is embedded in a medium with a dielectric constant e inside the simulation box. In addition, 
the box is filled with the neutralizing background charge. Periodic boundary conditions are 
applied. The corresponding electrostatic potential is determined from the Poisson equation 
with appropriate boundary conditions on the box boundary and ion surface.!! 



Figure 5. Finite-size dependence of the free energies of charging imidazole and imidazolium 
(filled circles and open squares on the right and left hand scale, respectively), as a function 
of the inverse volume of the simulation box, The top scale gives the number of water 

molecules. 
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Figure 6. Finite-size correction of the probability densities p{Uei) of the electrostatic 
energies Uei of Im(+). The uncorrected Uei histograms are shown with symbols, together with 
corresponding gaussian distributions. After correction for electrostatic and thermodynamic 
finite-size effects, the corresponding gaussian distributions "collapse" and agree closely for 
all system sizes of 16 < < 512 water molecules. 



Figure 7. Free energy of charging the polar imidazole Im(p) to the imidazolium cation 
Im(+) as a function of the order of the integration formula, (ijk) indicates the number i, j, 
and k of derivatives used at the uncharged, half-charged, and fully charged state.0 



Figure 8. Thermodynamic cycle with the free energies connecting the four states of the 
uncharged and charged imidazole and imidazolium as a function of a linear coupling pa- 
rameter. Shown are results for charging of imidazolium (solid line) , uncharging of imidazole 
(long dashed line), conversion of imidazole to imidazolium (dotted line) and conversion of 
uncharged imidazole to uncharged imidazolium (short dashed line). 



Figure 9. Comparison of the second-order perturbation expansion (dashed lines) with 
the reference free energies of charging Im(p) to Im(-|-) (top panel), uncharged Im(+) to the 
cationic Im(-|-) (middle panel), and uncharged Im(p) to the polar Im(p) (bottom panel). Also 
included as an inset in the bottom panel is a comparison with multistate gaussian models 
(symbols and dot dashed lines) shown with estimated statistical errors. The multistate 
expansions about the charged and uncharged states are shown with open squares and filled 
circles, respectively. 
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TABLES 



atom 


X 


y 


Q 


a 


e 


Imidazole 


Nl 


0.0000 


0.1105 


-0.090285 


0.325000 


0.71128 


C2 


-0.1091 


0.0282 


0.232373 


0.339967 


0.35982 


N3 


-0.0741 


-0.0983 


-0.715903 


0.325000 


0.71128 


C4 


0.0636 


-0.0984 


0.217356 


0.339967 


0.35982 


C5 


0.1120 


0.0298 


-0.374687 


0.339967 


0.35982 


HI 


-0.0009 


0.2112 


0.318027 


0.106908 


0.06569 


H2 


-0.2102 


0.0661 


0.102391 


0.242146 


0.06276 


H4 


0.1197 


-0.1905 


0.082346 


0.242146 


0.06276 


H5 


0.2119 


0.0700 


0.228383 


0.242146 


0.06276 


Imidazolium 


Nl 


0.0000 


0.1128 


-0.115106 


0.325000 


0.71128 


C2 


-0.1086 


0.0353 


0.010825 


0.339967 


0.35982 


N3 


-0.0663 


-0.0912 


-0.122786 


0.325000 


0.71128 


C4 


0.0719 


-0.0949 


-0.139642 


0.339967 


0.35982 


C5 


0.1140 


0.0344 


-0.122097 


0.339967 


0.35982 


HI 


-0.0018 


0.2141 


0.398875 


0.106908 


0.06569 


H2 


-0.2110 


0.0686 


0.230198 


0.242146 


0.06276 


H3 


-0.1274 


-0.1721 


0.402905 


0.106908 


0.06569 


H4 


0.1273 


-0.1872 


0.232002 


0.242146 


0.06276 


H5 


0.2131 


0.0766 


0.224826 


0.242146 


0.06276 


TABLE I: 


Coordinates x 


and y (in 


nm) and charges q 


(in elementary char^ 


;e units e) 



of the atoms in the planar imidazole and imidazolium from to the quantum mechanical 
calculations of Topol et al.il The Lennard- Jones parameters a and e (in nm and kJ/mol) 
are taken from the AMBER force field.!! Lorentz-Berthelot mixing rules were applied to 
combine the Lennard- Jones parameters of the solute atoms with those of SPC/E water.0 
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uncharged 


half-charged 


charged 


N 


ave 


var 


ave 


var 


ave 


var 


Imidazole 


16 


1.2 ±0.5 


209.6 ± 7 


-56.7 ±1.3 


436.8 ± 25 


-156.9 ±1.5 


495.6 ± 32 


32 


0.5 ±1.0 


224.4 ± 7 


-57. ± 1.6 


427.4 ± 22 


-158.7 ±1.5 


456.4 ± 26 


64 


1.8 ±0.6 


208.2 ± 7 


-58.8 ±2.0 


510.0 ±25 


-155.5 ± 1.5 


480.9 ± 15 


128 


1.0 ±0.7 


205.1 ±9 


-59.2 ± 1.6 


429.5 ± 16 


-156.1 ± 1.8 


491.9 ±20 


256 


1.3 ±0.5 


204.8 ± 7 


-59.2 ±1.6 


429.1 ±22 


-152.9 ±2.0 


473.4 ± 24 


512 


1.1 ±0.5 


209.6 ± 9 


-57.7 ±1.3 


410.8 ± 17 


-155.9 ±2.0 


539.2 ± 30 


Imidazolium 


16 


33.4 ±1.3 


1326.7 ± 10 


-226.5 ± 1.0 


1305.5 ±7 


-500.8 ± 1.0 


1391.0 ± 14 


32 


33.9 ±1.0 


1267.6 ± 12 


-215.0 ±0.9 


1267.0 ± 10 


-483.4 ± 2.0 


1375.0 ± 20 


64 


35.9 ±1.1 


1234.0 ± 12 


-208.7 ± 1.3 


1252.4 ± 13 


-473.9 ± 1.5 


1307.9 ± 18 


128 


34.9 ± 1.0 


1226.0 ± 14 


-203.5 ± 0.9 


1215.9 ±17 


-469.1 ± 1.3 


1291.2 ±22 


256 


34.1 ±1.1 


1222.5 ± 22 


-204.4 ± 1.2 


1256.3 ± 24 


-465.2 ± 1.5 


1330.2 ± 33 


512 


35.9 ± 1.5 


1248.8 ±18 


-201.2 ± 1.8 


1248.3 ± 26 


- 461.7 ± 1.5 


1305.3 ± 25 



TABLE II: Averages Ci (in kJ/mol) and variances C2 [in (kJ/mol)^] of the electrostatic 
energy of imidazole and imidazolium in the uncharged, half-charged, and fully charged states. 
Finite-size corrections have been applied. 
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simulation (kcal/e-mol) 




simulation (kcal/e-mol) 




-600 -500 -400 -300 -200 -100 

Uei (kJ/mol) 




2.5 3 3.5 4 4.5 5 
order of integration formula 



5.5 



